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Publication bias, the fact that studies identified for inclusion in a meta analysis do not represent 
all studies on the topic of interest, is commonly recognized as a threat to the validity of the results 
of a meta analysis. One way to explicitly model publication bias is via selection mo dels or weighted 
proba bility distributions. We adopt the nonparametric approach initially introduced bv lDear and Begg 
( 1993) but impose that the weight function w is monotonely non-increasing as a function of the p- 
value. Since in meta analysis one typically only has few studies or "observations", regularization of the 
estimation problem seems sensible. In addition, virtually all parametric weight functions proposed so 
far in the literature are in fact decreasing. We discuss how to estimate a decreasing weight function in 
the above model and illustrate the new methodology on two well-known examples. The new approach 
potentially offers more insight in the selection process than other methods and is more flexible than 
parametric approaches. Some basic properties of the log-likelihood function and computation of a p- 
value quantifying the evidence against the null hypothesis of a constant weight function are indicated. 
In addition, we provide an approximate selection bias adjusted profile likelihood confidence interval 
for the treatment efl'ect. The c orrespond i ng so ftware and the datasets used to illustrate it are provided 
as the R package selectMeta ( RufibacM 201l[) . This enables full reproducibility of the results in this 
paper. 
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1 Introduction 

Met a analysis has become a w idely used technique for synthesizing evidence from different studies, see 
e.g. Sutton and Higgin^ ( 20081 ) for an overview over recent developments. Publication bias, i.e. the fact 
that studies identified for inclusion in a meta analysis, do not represent all studies on the topic of interest, 
is commonly recognized as a threat to the validity of the results of a meta analysis . Overviews how 
to pre ven t, assess, and ad j ust fo r publication bias arc provided in Sutton et al. ( 2000l ). Macaskill et al 
(|200l[) . or lRothstein et all (|2005[) . 



Numerous tools to detect publication bias in meta analysis have been developed, see lRothstein et al.l (|2005 . 
Chapters 5-11) for an excellent overview of the current state-of-the-art. 

If one seeks to assess selection bias one typically requires some model for the sampling behavior of the 
observed effect sizes that expHcitly incorporates the selection process ( Hedges and VeveaL 2005 ). It is 
hence useful to distinguish two parts of such a model: the effect size part and the selection part. The 
former specifies what the distribution of the effect sizes would be if there were no selection whereas 
the latter explicitly models how the effect size distribution is modified by the sele ction process. Two 
differ ent classes of explicit selection models have been proposed so far for meta analysis ( Hedges and VeveaL 
20051 ). The first class depends on the effec t size estimat e , such as relative risk or odds ratio, and th e 
corresponding standard error separately, see Cooasl (ll999l ). iCoDas and Shil (l200nl) . ICopas and Shil (l200lh . 
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and the implementation in the R package copas ( Carpenter et al. . 20091) . This type of model is typically 
denoted "Copas selection model" . In the second class the weight is assumed to depend on the effect size 
only v i a, the jp- value ass o ciated with the study, see Hedged ( 19841 ) , llvenear and Greenhouse (ll988h . lHedg"esl 
(|l992l ). bear and Beg3 (|l992l) . iHedees and Veveal (120051 p. 149) or ICopas and MallevI (|2008l) for a test 
on selection bias that is robust against any form of selection function. The rationale to make the weight 
function depending on p-values, or equivalently on the standardized effect size, exclusively is that often, 
decisions about conclusiveness of medical rese arch results are based on statistical significance (only) . 
More specifically, following the development in lHedees and Vevea (l2005l) . let Y* be a random variable with 
density f{y\6, a) representing the effect estimate before selection, typically assumed to follow a normal 
distribution. Denoting the weight function by w(i/), the weighted density of the observed effect estimate 
Y is then given by 



f{y\0,a)w{y) 
j f{y\0,a)w{y) Ay' 



Whenever the weight function w is not constant, the sampling distribution of the observed effect size 
Y differs from that of the unselected effect size Y* and this difference, i.e. the shape of w, is a way of 
describing selection bias. 

Now, if larger values of Y* are more likely to be observed than smaller values, w{y) is a monotone non- 
decreasing function of the effect size y. Considering w on the scale of p-values this implies that w{p) as a 
function of the p- value is non-increasing, meaning that smaller p- values are more likely to be observed than 
larger p- values. In t his paper we propose a non- increasing estimate wij)) in the nonparametric normal 
model introduced by Dear and Beed ( 1992 ). Besides being a plausible assumption as elaborated above, 
nonparametrically estimating the weight function 'w{p) and imposing a monotonicity constraint has further 
advantages: 

• All parametric weight functions proposed in the literature are in fact non-increasing, see Section |3] 
for a brief discussion. 

• Typically, the number of studies that enter a meta analysis is small to moderate. For this reason, 
additional regularization, such as monotonicity, and therewith constraining the parameter space, 
may lead to more realistic but still fi exible esti mates of the weight fu nction compared to the purely 
nonparametric approach by [Hedges (1992) and Dear and Beed ( 1992 ). but without forcing a purely 
p arametri c mode l. This makes our approach less prone to misspecification. See also the comment 
in lHedged (119881 p. 118). 



• Restricting the parameter space, or shape of the function as in our case, typically yields estimates 
with better performance, e.g. measured in terms of me an squared error, if in fact the function to be 
estimated has the assumed shape, see e.g. iKellv (Il989l p. 937). 



• In con trary to e.g. kernel estimators or the penalized monotone estimator of iSun and Woodroofe 
(|l997| ) the estimator (u),0,(T^) defined below does not necessitate the choice of a smoothing or 
penalty tradeoff parameter (or a prior) and is therefore fully automatic. 

• Weight functions are primarily proposed as a n exploratory and infor mal means to assess the degree 
of publication bias which may be present, see Dear and Begd (1992, p. 240) or Sutton et al. ( 2000l 
p. 431). Specifically, if there is no selection effect at work, the former authors claim that the graphs 
of their estimated unconstrained weight function "provide visual confirmation of the lack of bias, 
demonstrating a seemingly random configuration of estimated weights." However, it is not without 
difficulty to identify th e model without biased selection from the estimated weight functions in sub- 
figures (a), (b), (c) of iDear and Beg3 ( 1992 . Figure 2). As reveal our examples in Section [71 the 
monotonicity assumption typically yields more insight in the actual selection process. 

In Section [2] we derive the log- likelihood function in our nonparametric model and provide some properties 
of it whereas in Section [3] we discuss different approaches to setup selection models and choose sensible 
selection functions. Section 2] elaborates on the computation of our proposed estimate. A discussion 
of statistical inference for the effect 6 and the random effects variance component are provided in 
Section [51 Specifically, in this section we sketch derivation of a profile likelihood confidence interval for 
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the selection bias adjusted treatment effect 9. A way to quantify evidence against the null hypothesis of no 
selection is described in Section [51 The paper is concl uded with the a nalysis of two well-known examples 
and a discussion of the software package selectMeta ( Rufibachl . l201l[) that enables full reproducibility of 
the results presented in this paper. 



2 The log-likelihood function and its properties 



To fix ideas, assume that there are n independent studies with normally distributed observed treatment 
effects Yi, i = 1, . . . ,n where ^(Vi) = 9 and Va.r{Yi) = rjf = uf + cP' . Here, itf is the known sampling 
variance in the i-th study (largely determined by the sample size in the i-th study and therefore considered 
known) and is a random effects component of variance representing the heterogeneity in the population. 
Typically, it is assumed that the effects follow a normal distribution, i.e. Yi ~ N{9,r]f) with realizations 
Hi. The two-sided p- values for the null hypothesis Hq : 9 = can then be compute d in each study as 
Pi = 2^{—\yi\/ui) and are, in accordance with the notation of Dear and Beed ( 1992[) . considered to be 
ordered and denoted by p„, . . . where p„ is the smallest and pi the largest. Furthermore, let Pn+i = 
and po ~ I. Assume that the selection process is governed by the non- negative weight function w that 
assigns to an effect estimate the likelihood that it is observed. Then, the likelihood function of the observed 
effect sizes y = (yi, . . . , y„), given the weight function w, the quantities 9, ct^, and u = . . . , it„), 
amounts to 



L{y\w,9,a^ ,u) ~ P(?/i|i-th study is published) 



i=l 



where we introduced the normalizing constant 

Ai{w,9,a^,u) = 



Ai{w,9,a'^,u) 



(1) 



(2) 



and $ as well as (p, the cumulative distribution and density function of a standard normal distribution. 
Now observe that in ([1]) the unknowns are the function w and the parameters 9 and . Numerous 
suggestions have been made to estimate these unknowns, where these proposals differ by the assumptions 
they impose on the selection function w, s ee th e discuss i on in Section [31 



Dear and Beg3 (Il992l ) and 
left-continuous step function of the p-value. In 



In this paper, as in 



Hedges 




we posit that the weight function if is a 
the discontinuities of w are fixed at, say, 



iHedge. , _ 

"psychologically motivated" values, whereas iDear and Begd (1992) group the p- values in pairs and assume 
equal values of w for two adjacent observed p- values. Here, we adopt the latter approach noting that the 
former model fits in our framework equally well. More specifically, the weight function is, for p e [0, 1], 
defined as 



'){p) 



if P2J-2 >P> P2j 

J Pn-1 > P > if 71 odd 
I Pji > P > if rt even 



where j = 1 + L'/2J = 1, . • . and k is the number of categories that are built from the initial p- 
values through pairing. For reasons of identifiability one is not able to set up a li kelihood assuming a 



piecewise constan t weight function with out some sort of grouping of the p- values (see Sutton et al. , 2000l . 



Section 2.3.7 and Dear and Begd . 19921 . Section 2). Fo r the description of a "pur e selection model" and 
the necessary modifications of the problem, we refer to ISun and Woodroofd (Il997h . 
The weight function on the scale of the outcomes y writes as: 



w{y) 



Wjl{-u,<p-\p2j/2) > \y\ > -u,$-i(p2j_2/2)} for i = 1, 



(3) 



To see this, note that if the p-value in study i gets the weight Wj assigned, this p-value is computed for 
a test statistic \yi\/ui and equal to ph = 2$(— |yi|/ui), what gives \yi\ = —Ui^^^{pii/2). Plugging in 
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this form for the weight function w into ^ and taking the fog yields the final weighted log-likelihood 
fu nction l(w,9,a^) for t he parameter vector {w,9,a'^) € R'"""''^. This log-likelihood was initially derived 



Dear and Beed (1992). However, here and in the appendix we summarize its detailed development and 



discuss some properties and additional computational facts. For the log-likelihood function we get 

k 

A ; lnp'7/i, — N ^ focr n, — ^ 



k n A ^ fi 2 

l{w,9,a') = _(„/2)log(2^) + ^A,logz.,-^logr;,--^(^i— ) -J^logA 



where Ai,i = 1, . . . , n are the normalizing constants defined in ([2]). Straightforward computations for any 
c > yield that the log-likelihood function can be written as 



l{cw,9,cr^) = l{w,9,cr^)+\og{c)\^(j2^j) 



3 = 1 



= l{w,9,a^) + \og{c){Xi^l). 



(4) 



Two important observations can be made for I. First, the quantities Xj should, in principle, correspond 
to the number of p-values in any interval (p2j,P2j-2], j ~ ^, ■ ■ ■ ■, k, i.e. Ai = 1, Xj = 2 for j = 2, . . . , fc — 1 
and Xk = 1 -I- 1 {if n is odd}. However, the choi ce Ai = 1 would imply by (|1]) that the maximizer w of / 
was not identifiable. To overcome this problem. Dear and Beed (1992, P- 239) advise setting Ai = 2, so 
that ([4]) simplifies to l{cw, 9,(t^) ~ l{w,9,a'^) + logc, making w (1) identifiable but (2) estimated with a 
slight negative bias. For reasons of simplicity we choose Ai = 2 in the examples in Section [71 In the code 
collected in selectMeta (jRufibachl . 120111 ) the weight Ai can be set to an arbitrary value. 
Second, ^ entails that we must have w <£ (0, 1] since lii = 1 once Ai > 1 is chosen. To see this, assume 
w with wi < 1 the largest element and choose c = l/wi > 1. This yields 



l{cw,9,a^) 



l{{l,W2/wi, . . .,Wk/wiy,9,a'^) 
l{w,9,a^)-log{wi){Xi-l) via® 



and thus 



l{{l,W2/wi,...,Wk/wiy,9,a^) > l{w,9,a'^) 



if Ai > 1. However, as discussed in Dear and Beed ( 19921) . the actual selection probability is typically less 
than 1 for all studies under consideration since some selection is going on for all studies, or rather the 
corresponding p-value. As a consequence, the estimated weights are only relative. Since no information 
is available on the p-values of the unpublished studies, one is not able to estimate the weight function 
directly. 



The primary goal of this work is to specialize the approach of iDear and Beed (jl992f ) to a monotone 
selection function w. Thus, we followed the framework developed in the latter paper, to enable straight- 
forward comparison of the newly introduced monotone weight function to the existing approaches and 
thus made the weight function w(p) depending on two-sided p- values. However, the entire framework can 
straightforwardly be adapted to one-sided p-values. 

Ideally, in order to apply standard algorithms to maximize a log-likelihood function one appreciated if it 
would be nicely behaved, i.e. strictly concave and coercive. Unfortunately, this is in general not the case 
for l{w,9,a^). Instead, plots of Z as a function of one of its arguments reveal that it is not necessarily 
concave in Wj,j = 1, . . . ,k and a. However, these same plots strongly indicate that / is at least unimodal 
with a unique maximum, although we are not able to provide a formal proof of this property or some 
(even stronger) surrogate, like e.g. log-concavity of I. Assuming that in fact / were unimodal. Lemma [TT] 
below would then imply that a maximizer always exists. In addition, the expression of the likelihood in 
the lemma also sheds some light on the peculiar structure oi I. To state Lemma ETTl let p = {w,9,(j^). 

Lemma 2.1. Assume that n > 3, Wk > 0, and Xj < n for all j. Then, the log-likelihood function I is 
continuous as a function of {w,9,a^) and coercive when one or more coordinates approach the boundary 
of the domain, i.e. if \\p\\ — > oo and/or if Wj — > for at least one j < k, then I — )■ — cxd. 

Note that / remains finite if cr and all other arguments are kept fixed. 
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Proof of Lemma 12.11 First, note that for a fixed i, 

log Ai = \og{wiH.a + . . . + WkHik) 

= log(ti!i) + log(i7ji + . . . + WkHik/wi) 

i=i j=i 
Using this, the log-Ukehhood function can be written as 

k n 1 ^ 2^ 

l{w,0,a') = _(n/2)fog(2^)+^A,log«;,-^log77,--^(^^) -^logA, 

k 1 ^ 

= -in/2) log(2^) + ^(A, - n) \ogw, - - ^ \og{u^ + a^) 

^ n n k _^ 

-^Y.(y^-enu^ + aT'-Y.\og{Y^H,,{ul,,^^wiy }. (5) 

i=l i=l j=l 

Let be a sequence of vectors such that \\pj,\\ — > cxd as r — > oo. From the definition of Hij in Appendix El 
it is clear that Hij G [0, 1]. The assumption Wk > entails that at least one Hij is different from 0. From 
(O it is then not difficult to see that 

l{wr,Or, {cr'^)r) " > ~oo as r ^ oo 

for either combination of possibilities, i.e. Wrj oo for one or more j's and/or \6r\ — ^ oo and/or 
(cr^)r- — >■ oo. Representation ([5]) also illustrates the continuity of I. Now, from ([5]) we can derive that 

k 1 ^ 

l{w,9,(j'^) = -(n/2)log(27r) + Ailogw;i+^(Aj- -n) log - -^log(u.? + cr^) 

i=l 1=1 j=2 

Without loss of generality assume that ||p,,|| — > oo or to some constant, but that Wr,i 0. The above 
representation then readily implies that l[wr,9r, — oo. □ 



3 Monotone selection function 



Fo r a thorough review of weight functions w proposed in the literature for meta analysis we refer 



to 



Sutton et al.l (120001 Sectio n 2). The spectrum ranges from (1) fully par ametric proposals as in 



vengar and Greenhouse! ( 1988[ ) , the weight fu n ction proposed in the comment to Iyengar and Greenhouse 
19881 ) by Hedges, or tho se in Preston et al.l (12004 Section 3.2) to (2) nonparametric models as those 



discussed in lHedges ( 1992 ) and Dear and Begd (Il992). Many of th ese functions have also been co nsidered 
in a Bayesian framework, see the discussion in ISutton et al. I ()2000l . Section 2) or lSillimanI (Il997allbl) . 
In general, there is little empirical evidence to guide the choice of weight functions (Hedges, 19881 p. 119). 
However, the literature generally agrees that weight functions that depend on t he effect size through the 
corre sponding p-value should b e non-incr easing as a function of the p-y a.lue, see Dear and Begg ( 19921 p. 



2381. Ilvengar and Zhaol (|l994L p. 38), or II^ (|200l[) . In the review bv ISutton et all (l2000l Section 2.3) 
all eigh t weig ht functions that are discussed are in fact monotone non-increasing. On the other hand, 
[Hedged ( 1992 . p. 249) argues that "It is probably unreasonable to assume that much is known about the 
functional form of the weight fu n ction. " Combining these two demands we therefore propose to adopt 
the approach bv iDear and Begg ( 1992h . i.e. to use the log-likelihood function l{w,9,a'^) developed in 
Section [21 but maximize it over the set 



V 



{{w,e,a'-) : 1 



Wk 



> 



,<7>0}. 
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So, we aim at computing 



{w,9,a'^) = argmax l{w,6,a'^). 
For completeness, we also state the unconstrained problem of iDear and Beed (|l992l ) 



argmax l{w,9,a ). 

{wM,a-^)e[0,l]'' xRx [0,00) 



(6) 



(7) 



To conclude this section we would like to point out Givens et al.l ( 1997L p. 228), for two reasons: First, to 
the best of our knowledge these are the only authors who explicitly estimate a monotone non-increasing 
weight function. However, in a Bayesian context via rejection sampling. Second, they remark that "Such 
a [monotonicity] constraint is much harder to put in place in the frequentist setting...". Here, we close 
this gap for the likelihood setup described above and provide corresponding R software, see Section [51 

Sensitivity virith regard to assumptions. As discussed in the introduction, explicitly modeling the 
selection function depending on the study p-value is only one way of adjusting for selection bias, the 
most prominent alternative being the Copas selection model that ma kes sele c tion depending on the effect 
size estimate T and t he corresponding s t andard error a separatelv (ICopa 1 11999I ICopas and Shil . I2OOOI 



. ^ . pJCo 

Copas and Shi l200l[ [Hedges and Veveal . l2005i ICarpenter et al.l . 120091 ). Howe ver, making the selection 



function depending on the p- value only has the longest history in meta analysis (jHedges and Veveall2005 . 
p. 149). A potential constraint of the latter models is that they treat equally significant results in the 
same way, irrespective of the size of the underlying study and the direction of the effect. Thus, a potential 
next step to generalize the approach proposed here is to setup a two-dimensional selection function that 



• non-increasing as a function of p- values and 

• non-decreasing as a function of the underlying study size. 

A potential source of misspecification is an inappropriate choice of w's shape. However, as elaborated 
in Section [TJ monotonicity seems a very plausible assumption for a selection function and all proposed 
parametric approaches are in fact non-increasing. Finally, in order to correct for publication bias, selection 
models must substitute assumptions for data that are missing. In our scenario, we stipulate that the 



form of the unselected distribution of the effect size estimates is normal. However, Hedges and Vevea 



(jl996f ) performed a large simulation study assessing robustness of estimation from selection models to 
misspecifica tion of effect distribution and concluded that, surprisingly, the procedure is rather robust in 
this regard (jHedges and Veveal . [2OO5I) . 

Note that estimation in the C opas selection mode l is not free of difficulties and estimation can be impossible 
for certain parameter values ([Hedges and Veveal [20051) . 



4 Computational aspects 



Having formulated the problem ^ it remains to numerically compute {w,9,&'^). As a consequence of 
the considerations in Section [21 properly maxi mizing I is somewhat delic ate, even when looking at the 
unconstrained problem ([7]). Note that neith er Dear and Begg (.199^ nor Hedged ( 1992 ) discussed this 
aspect. They both mention ( Dear and Beggl 19921 p. 240 and Hedges . 19921 p. 251) that a multivariate 



Newton- Raphson procedure can be used to find the unconstrained maximu m . . . ,vjk*,0*,o''^)- To 

avoid inversion of the cor respond i ng He ssian matrix [Dear and Begd ([1992I) use an EM-type algorithm 
which is also discussed in Hedged ( 1992 ). Namely, one iterates optimization for one entry of {w,0,a'^) 
at a time using Newton-Raphson until convergence. We implemented both approaches and surprisingly, 
the one-entry-at-a-time version turned out to be more stable and quick enough and has therefore been 
implemented in selectMeta, see Section [51 for details. 

However, we did not find a way to generalize this approach to find our new, constrained estimator (iD, 9, ct^) 
defined in ([6]). In fact, to find this estimator we have to maximize a (most likely) unimodal, coercive but 
generally non-concave function under constraints, a non-trivial global optimization problem. The solution 
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wc present makes use of the so- called evolutiona r y glob al optimization via the differential evolution (DE) 
algorithm, initially described in lStorn and Pried ([1997(). This algorithm is particularly well-suited to find 
the global optimum of a real- valued function of real- valued parameters, such as our log-likelihood function 
I. Neither continuity nor differentiability is a necessary property of the target function that is maximized 
by a DE algorithm (as a m atter of fact, our I above is di f ferent iable) . An implem entation of a DE algorithm 
is made available in the R ( R Development Core Team . 2010t ) package DEoptim ( Ardia and Mullen . 201Clt ). 
The function DEoptim allows for unconstrained global maximization. To account for constraints, as in our 
case given by the monotonicity assumption in the parameter set "P, one simply integrates the constraint 
within the function to optimize by penalizing deviations from the constraints with — oo. Set up this way, 
the function DEoptim quickly and reliably delivers the maximum of I over V. For a description of the 
implemented software we refer to Section [S] 



5 Statistical inference on 9 and cr^ 



We agree with lHedeed (jlQSSL p. 120) when he says that "Although I am enthusiastic about the development 
of varied and realistic models for estimation under selection, I do not believe that estimates from any one 
of these models should be taken too seriously." This approach of considering selection models a way of 
ex ploring the selection m echanism, but not to estimate the parameters 9 and ct^, is further supported 
by iDear and Begg ( 19921 . p. 240) claiming that "The procedure presented here [meant is their selection 
model] is intended primarily as a means of informally exploring the degree of publication bias which may 
have operated in the selection of studies contri buting to a meta ana lysis. Inference about 6 and should 
be considered secondary at this stage." or by ISutton et al. I (I2OOOI p. 431/439): "Clearly, it is far more 
desirable to alleviate the problem of publication bias rather than try to model it analytically. [...] Hence, 
the weight function obtained is used to provide a visual display of the relative weight function for the 
purposes of identifying publication bias, and is not used to adjust the pooled estimate." If selection bias is 
suspected from looking at w (or w^) one should primarily focus attention on the possible causes of bias, 
e.g. initiate a search for "missing" studies, rather than using the model to adjust 6 and (t^ for publication 
bias. 

However, to get a complete picture we would like to ske tch a way of making select i on bi as adjusted 
inference for 0. As elaborated in the seminal paper by iMurphv and van der Vaa"rt (I2OOOI) on profile 
likelihood in presence of an infinite-dimensional nuisance parameter, ordinary profile likelihood inference 
may still be applicable if the entropy of the function class of the nuisance parameter is not too large. 
When seeking inference for 9, the nuisance parameters are a and the estimated monotone weight function, 
w. That the cla ss of monotone functions is not "t oo large" in ter ms of entropy and th us access i ble for 
the approach by iMurphv and van der Vaart (j2000h is discussed in iFan and Wonel (120001) . IChosliI (|200it 
uses this approach to provide inference on a one-dimensional parameter with an estimated monotone 
nuisance func tion in the evaluation of a biomarker. Here, by appealing to the above profile likelihood 
arguments of Murphy and van der Vaart ( 2000| ). we get that the likelihood ratio-based statistic for the 
parametric component will have a limiting distribution with one degree of freedom. Based on this 
result we derive a selection bias adjusted profile likelihood confidence interval for 9 which is implemented 
as the function DearBeggMonotoneCItheta in selectMeta. Note that this procedure is approximate and 
rigorous theoretical justification of this approach will be provided elsewhere. 



6 Quantifying the evidence against a constant weight function 

Obviously, one would like to have a mean to quantify the evidence against the null hypothesis of a 
constant weight function w. The only reference we are aware of that deals with a similar problem is 
IWoodroofe and Sun (I1999I) . A monotone if is considered, but the density of the effect sizes is assumed to 
be entirely known what precludes application to our situation. 

However, as an alternative we suggest a simulation procedure to get a p- value that quantifies the evidence 
against a constant weight function w, based on our new monotone estimator. Before describing compu- 
tation of such a p- value let us introduce the density function g of the distribution of p- values for a meta 
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analysis with true effect 9, variance u^, and random effect component a^: 



{-a$-i(p/2)-^}/r/ {a<^-Hp/2)-e}/rj 



0{<J>-ib/2)} 



(8) 



where rf ^ v?' + . This is tlie density generated by a test o f the hypothesis JJp : Y ^ A^(0,(T^) vs 



Hx:Y ^ N{e,rf). Note that / simplifies to denoted by g[p) in iDear and Beed (|l992l p. 240) for their 
choice u = 1, (T = of parameters. 

The log-likelihood function I does not depend on the p-values only, but also on the sign of the initial 
effect size y. So when simulating p- values from /, to be able to compute the log-likelihood function, 
it is therefore not sufficient to generate a sample of p-values from the density / only (via numerical 
inversion of the quantile function corresponding to /) but one also needs to randomly generate the signs 
of the corresponding "observations" y. For each generated p-value p and a Bernoulli random variable 
Z ~ Ber(l/2) we therefore compute y* — — 7i$^^(p/2) and set 



y 



y* if Z = 0, 

26 -y* if Z = 1. 



Note that to simulate a p-value from a distribution with density one could equivalently first generate 
a random number y drawn from the distribution N{9,rf') under the above alternative and then compute 
the p-value as p = 2<i>(— |j/|/it). 

Now, to generate a (one-sided) p-value for the null hypothesis of a constant weight function we proceed 
as follows: 

1. As test statistic to assess constancy of a monotone weight function w we choose T ~ minw. 

2. Compute estimates 9^ and (Tq from the observed collection of p- values pi,...,pn in a standard 
random effects model. Also compute the monotone weight function wo based on this collection. 

3. Draw samples (pji, . . . ,Pjn) of p- values for j = 1, . . . , M where pji follows a distribution with density 
f{-;9o,a-o, \/uf + d^) for i = 1, . . . , n. For each of these samples also compute the monotone weight 
function ilij. It is important to realize that these samples of p- values come, by construction, from 
the null model, i.e assuming no selection bias. 

4. Compute the test statistics Tq = minwo and f — min Wj for j — 1, . . . ^ M. 

5. The proposed approximate p- value that quantifies the evidence against a constant weight function 
is then 

1 + #{J<A/ : fo<rW} 

P = TTTi ■ 



The function DearBeggMonotonePvalSelection implements this procedure in selectMeta. 



7 Examples 



Op 

ofll 



en classroom data 



Dear and Begg 



As a first example, and to compare the monotone to the non-monotone approach 
we re-analyze the famous open classroom education data i nitial ly presented by 



Hedges and Olkin (|l985l p. 303) and re-analyzed bv llvengar and Greenhous3 (1988) and Dear and Begg 
(|l992l) . For convenience, the data is reproduced in Table fTl (compare llvengar and GreenhouseLll988[ Table 
4). 

All these studies assessed the effect of open vs. traditiona l education on student creativ ity, measured by 

some continuous quantity (in fact, we did not find neither in lHedges^ and Olkin . 19851 nor in llvengar and Greenhous^ . 
19881 the exact description of what was actually measured) . In Table [2 Ni denotes the sample size in 
each of the two samples (so all the studies were perfectly balanced), yi is the mean difference (the effect 
measure), are the standard errors, and pi are the computed p-values, pi = 2^(—\yi\/ui). 
In Figure [T] we present the following estimates of the weight function: 
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i 


Ni 


Vi 


Ui 


Pi 


1 


10 


0.081 


0.45 


0.86 


2 


10 


0.308 


0.45 


0.49 


3 


39 


-0.178 


0.23 


0.43 


4 


50 


-0.234 


0.20 


0.24 


5 


10 


0.598 


0.45 


0.18 


6 


22 


0.563 


0.30 


0.06 


7 


40 


0.535 


0.22 


0.02 


8 


36 


0.779 


0.24 


0.0009 


9 


20 


1.052 


0.32 


0.0009 


10 


90 


-0.583 


0.15 


0.0001 



Tabic 1: Studies of effects of open vs. traditional education on creativity. 



• The parametric weight functions wi and W2 proposed in Ivengar and Greenhousg ( 1988 . p. 113), 

• the nonparametric variant of iDear and Beed ( 19921 ) . 



the nonparametric variant of iDear and Beed (|l992r ). 
• and our new proposal: the nonparametric weight function constrained to be non-increasing. 



As in lDear and Bead ()1992i) we provide two plots: one with the original p- value scaling of the x-axis and 
one where on the x-axis we plot the limits of the pairwise groups of p- values, where these limits are equally 
spaced. Note that in the latter plot (1) the parametric weight functions are not displayable and (2) one 
must carefully study the horizontal axis to determine the p-values represented by the estimated weight 
functions. At the bottom of the first plot, we also indicated the observed 10 p-values with vertical ticks. 
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Figure 1: Estimated weight functions for the open classroom education dataset. 



Since the nonparametric estimate of iDear and Bead (|1992I ) is already "almost" non- increasing it comes 
without surprise that the monotone estimate of w is very similar to its unconstrained co unterpart. The 
estimated weight functions clearly indicate publication bias, an observation already made in lDear and Begg 
(|l992l . p. 243). This IS further supported by the p- value computed according to the procedure outlined in 
Section |6] that amounts to ^education = 0.096 (for M = 1000 runs). Now, having an estimate of w at hand 
yields some more insight in the selection process: According to the monotone estimate, the probability of 
a p-value that is larger than 0.001 to be published is only 28.0% compared to a p- value < 0.001. 
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Finally, estimates for this datasct from our monotone selection model are 9 = 0.14 and ir^ — 0.11 with 
95% approximate profile likelihood confidence interval for 9 of [—0.08,0.57]. Compare these to estimates 
received from a standard random effects model that amount to 9 — 0.26 and = 0.30 with estimated 
confidence interval for 9 of [—0.12,0.65]. The adjustment for selection thus attenuates the effect estimate 
and narrows the confidence interval but does not change the conclusion about significance of at a 
significance level of a = 0.05. 



Envi ronmental tobacco smoke data In our second example we discuss a mcta analysis (jHackshaw et al 



19971 ) of 37 studies concerned with the effect of environmental tobacco smoke on lung-cancer in lifetime 
non-smokers. The effect in these studies is quantified via the log relative ris k. Whether this meta analysis 
suffers from publication bias has been a matter of ongoing controversy, see iRothstein et al.l (j2005l p. 91). 
In the original publication a peculiar form of "failsafe iV" analysis was conducted and the autho r s con- 
cluded that there is no reason to suspect publication bias. In a re-analysis however, Cooas and Shi (j2000[) 
(see al s o the correspondence following that paper on www.bmj .com), applying the method introduced in 
Copad (jl999h . came to the conclusion that "the possibility of publication bias cannot be ruled out alto- 
gether, and at least some publi cation bias is nee ded to explain the t rend w e found." However, neither the 
funnel plot, nor the method by[Co2ai([i993), or lCopas and MallevI (l2008l) yields real insight in the nature 
of the selection process that may be at work. On the other hand, we can estimate the weight function via 
the unconstrained and t he monotone approach , see Figure [51 

First, unlike claimed in Dear and Begd ( 1992 . comment to Figure 2), note that from the unconstrained 
estimate it is not evident whether publication bias is operating on this dataset. 

Again, we can gain some insight in the possible selection mechanism by looking at the estimated weight 
function which reveals an interesting pattern: One observes four distinct regions which are given by the 
intervals [0,0.03], (0.03,0.17], (0.17,0.77], (0.77,1.00] where w is constant. These regions are indicated 
with vertical dashed lines in the left plot of Figure [2] Not surprisingly, sharp drops in the weight function 
appear around "psychological barriers" for p-values, namely 0.05 and maybe 0.15. In passing we remark 
that discontinuities of the estimated weight function can only happen at actually observed p- values. 
The probability of selecting a study with p-value larger than 0.17 is only 64.8% of that of one with a 
p-value at most 0.17. In addition, the relatively small p-value of Ptobacco — 0.13 computed according the 
method described in Section [5] reveals some evidence against a constant weight function in Figure [2] For 
these reasons it seems th erefore plausibl e that publication bias is at work here and we thus agree with the 
conclusion of ICopas and Shi ( 20001) and [Hedges and Veveal (j2005l p. 164). 

Furthermore, in a standard random effects meta analysis model, we get estimates 9 ~ 0.21 and = 0.02 
with estimated confidence interval for 9 of [0.12,0.31]. These estimates are attenuated to 9 — 0.17 and 
(T^ = 0.01 in the monotone selection model, with 95% approximate profi le likelihood confidence in terval 
for 9 of [0.08,0.26]. These changes are very similar to those observed bv [Hedges and Vevea (2005) when 
choosing their somewhat related weight function. And the significant effect of environmental tobacco 
smoke on lung cancer persists after adjusting for selection bias. 

Finally, to illustrate the computation of the suggested p-value, in the lower plot in Figure [2l we have 
plotted in grey the estimated weight functions for the M = 1000 samples generated under the assumption 
of no selection. This gives an impression what selection function can be expected under no selection. 



8 Software and reproducibility 

Although some of the methods discussed here have been around for some tim e, it seems as if none of them 
has found its way in the daily routine of meta analysts. ISutton et al. I (|2000l p. 439) expla in this lack of 



use by "One of the explanations for this is almost certainly the complexity of many of the approaches 
(particularly the selection models, and Copas' approach), and the lack of user- friendly software available 
to implement them." In addition, "A further reason for low penetration is possibly lack of acceptance of 
such methods, and "Selection models are quite sophisticated and there is currently a lack of software to 
implement them" . To foster the use of selection models by meta analysts we have therefore implemented 



the parametric weight functions wi and 1^2 from llvengar and Greenhouse! (|l988f ) as well as maximum 
likelihood estimation of their corresponding parameters, 

the nonparametric method of lDear and Beggl(|l992l) . 
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Figure 2: Estimated weight functions for the environmental tobacco smoke dataset {n = 37, fc = 19). 
Thick Hues: monotone non-increasing weight functions. Thin hues: Unrestricted weight functions. Dashed 
vertical lines at 0.03, 0.17, 0.77. 



• our new variant that provides a monotone version of the latter estimate, including estimation of the 
selection bias adjusted estimates of 9 and cr as well as the approximate profile likelihood confidence 
interval for 9, 

• the density, distribution, and quantile function as well as random number generation from the p- value 
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density ©, 



• and the procedure to compute a p-value to assess the null hypothesis of no selection introduced in 
Section [HI 



in a new R package selectMeta (jRufibach, which is available from CRAN. In addition, we provide in 

selectMeta the two datasets analyzed in Section [71 

Making the software and datasets discussed in this paper accessible enables reproducibility of the results 
and plots. The code that generates Figures [1] and [2] as well as the computation of the p- value introduced 
in Section[n]for these two examples can be f ound in the h elp file fo r the function DearBegg in se lectMeta. 
This document was created using Sweave (jLeischl . [20o3) . (jKnuthl 11984 iLainportl . Il994[) . R 2.12.2 

( R Development Core Team . 2010f) with the R packages selectMe ta ( RufibachT 20 1 ll Version 1.0.3), 
DEoptim (lArdia and Mullenl . I2OIOI Version 2.0-9), met a dSchwarzeil I2OIOI Version 1.6-1), reporttools 
(jRufibachl . I2OO9I Version 1.0.5). and cacheSweave (jPentd . l2008l Version 0.4-5). 



9 Final remarks 

We propose and analyze a new type of monotone frequentist nonparametric weight functions as a visual 
tool to gain insight in the study selection process when publication bias in meta analysis must be suspected. 
Selection models have not yet entered the standard toolbox of meta analysts, presumably due to lack of 
easy accessible software. Our goal was to redu ce this gap by co llecting many existing and our new approach 
as functions in a new R package selectMeta (jRufibachl . 120111) . 

More research is necessary to popularize selection models. We intend to develop a smooth version of our 
new estimator by imposing not only a monotonicity but also a smoothness constraint on the log-likelihood. 
However, already difficult algorithmic aspects arc not facil itated by such additi o nal re gularization struc- 
ture. We further plan to apply and adapt the method of ISun and Woodroofd (|l997l ) to meta analysis. 
Finallv. [Hedges and Vevea ( 2005 . Eq. 9.6) describe how to incorporate not only as a simple number, but 
rather depending on covariates in a regression model. This approach should also allow for a generalization 
to other than the specific selection model they look at. 
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A Derivation of relevant quantities 

In this brief appendix, we provide some additional computations that lead to the log-likelihood function, 
merely for the reader's convenien ce. As a matter of fact and since the log- likelihood used in this paper is 
the one of Dear and Beggl(|l992l) . a derivation of I can also be found there. However, here we try to be a 
bit more explicit. 

Since the weight function w is defined to be piecewise constant with values Wj , the normalizing constants 
Ai simplify to 



A, 



E 



i=i 

k 



y:Wi(y)=Wj 



0(^) 
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The quantities Hij can be computed as follows: 



0(^) d. 



dy 



bi.2j-2<\y\<bi.2j ^ 



where we defined 



$(ay)-$(fey) + 3'(c.j)-*(t^u) 



_ Ui\y2j\/U2j - 9 , _ Ui\y2j-2\/u2j-2 - 9 



see 



_ -Ui\y2j-2\/u2j-2 - , _ -Ui\y2j\/U2j - 9 

Cij dij - , 

■Hi Vi 

Dear and Begg ( 19921 Appendix). Consider the following "boundary cases": Defining po = 1 and 



P2k = 0, we get bi^o = and bi^2k = oo, what immediately entails 



Ui\y2\/u2-9\ -u.,\y2\/u2 



^ ^^ u.\y2\/U2-u ^ _^^ -u.,\y2\/U2-t^ ^ ^ $(a,) - <i>Ki). 
On the other hand, 



^ 1 _ ^(U.\y2k-2\/u2,-2 9^ ^^^-U,\y2k-2\/u2k-2 9^ ^ + 

Plugging all the above quantities into ([T]), the log- likelihood function amounts to 
l{w,9,a^) = logL{w,9,a^) 

n n ^ n k 

= Y: iog^.(y.) + E iog{,r v(^) } - E i°g(E -.^».) 

= -|log(27r) + EAjlog^«,-Elog''^- 2^(^^) -El°g^- 

j=l i=l i=l 1=1 
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